
#include <stdlib.h>
#include <stdio.h>
#include <gmp.h>
#include <time.h>



void gauss_legendre(unsigned precisionBits, mpf_t pi, unsigned int *i);

int main(int argc, char *argv[])
{

	unsigned prec_bits = 35000000;
	time_t start, end;
	double dif = 0.0;
	unsigned int i;
	mpf_t pi;
	mpf_set_default_prec(prec_bits);
	mpf_init(pi);

	time(&start);
	gauss_legendre(prec_bits,pi,&i);
	time(&end);

	dif = difftime(end,start);

	printf("Time: %2lf\n",dif);
	printf("Iterations: %u\n",i);
    printf("PI = \n");
    mpf_out_str(stdout, 10, 10000001, pi);
    printf("\n");
    mpf_clear(pi);
	return 0;
}

void gauss_legendre(unsigned precisionBits, mpf_t pi, unsigned int *i)
{
	unsigned p;
	mpf_t a;	
	mpf_t a_aux;
	mpf_t b;
	mpf_t b_aux;
	mpf_t t;
	mpf_t temp;

	/*================== Init variables ==================*/
	mpf_init(a_aux);
	mpf_init(b_aux);
	mpf_init(temp);
	mpf_init(a);
	mpf_init(b);
	mpf_init(t);
	
	
	/*================ Initial Conditions ================*/		
	mpf_set_ui(a,1); /* a = 1 */

	
	mpf_sqrt_ui(b,2); /* b = 1/sqrt(2) */
	mpf_ui_div(b,1,b);

	
	mpf_set_ui(t,4); /* t = 1/4 */
	mpf_ui_div(t,1,t);

	p = 1;

	for(*i = 0; *i < 22; (*i)++)
	{

		/*A(n+1)*/
		mpf_add(a_aux,a,b);
		mpf_div_ui(a_aux,a_aux,2);

		/*B(n+1)*/
		mpf_mul(b_aux, a, b);
		mpf_sqrt(b_aux, b_aux);

		/*T(n+1)*/
		mpf_sub(temp, a, a_aux);
		mpf_pow_ui(temp, temp, 2);
		mpf_mul_ui(temp,temp, p);
		mpf_sub(t,t, temp);

		mpf_set(a, a_aux);
		mpf_set(b, b_aux);
		
		/*P(n+1)*/
		p <<= 1;

		mpf_add(pi,a,b);
		mpf_pow_ui(pi,pi,2);
		mpf_mul_ui(temp,t,4);
		mpf_div(pi,pi,temp);
	}

	mpf_clear(a);
    mpf_clear(b);
    mpf_clear(t);
    mpf_clear(a_aux);
    mpf_clear(b_aux);
    mpf_clear(temp);
	
}